library(xtable)
library(foreign)
library(countrycode)
library(Design)
library(apsrtable)
library(MASS)
library(car)
library(pwt)
library(lattice)
library(effects)
library(ggplot2)
#library(psychometric)


data<-read.csv("data_for_export.csv",header=T,as.is=TRUE)

## OPTIONALLY, GRAPH SOME THINGS
## source("graph.r",echo=TRUE)


#############################
##
## Check some correlations
##
#############################


cor(data$tor,data$vul,use="pairwise")
cor(data$defacto_t,data$magetti,use="pairwise")
plot(data$defacto_t,data$magetti,xlim=range(data$defacto_t,na.rm=T),ylim=range(data$magetti,na.rm=T),type="n")
abline(lm(magetti~defacto_t,data=data))
text(data$defacto_t,data$magetti,paste(data$iso3cc,data$sector,sep="-"),cex=0.4)

my.steps<-7
working<-data[!is.na(data$Ourscore),]
working$fm<-as.numeric(working$sector=="Financial markets")
working$energy<-as.numeric(working$sector=="Energy")
working$Laver.Hunt[is.na(working$Laver.Hunt)]<-mean(working$Laver.Hunt,na.rm=T)
lm.mod1<-lm(defacto_t~Ourscore+rol+runvetops+runpolarization+log(runpop)+log(rungdpch)+log(age)+HallGingerichCoordination+Laver.Hunt+Turnover+Staff,data=working) #,x=T,y=T)

lm.mod1$se<-robcov(ols(formula(lm.mod1),data=working,x=T,y=T),working$iso3cc)$var


lm.mod2<-step(lm.mod1,steps=my.steps)
#lm.mod2<-update(lm.mod2,.~.-age)
lm.mod2$se<-robcov(ols(formula(lm.mod2),data=working,x=T,y=T),working$iso3cc)$var

## NOW WITH TOR
lm.mod3<-update(lm.mod1,tor~.)
lm.mod3$se<-robcov(ols(formula(lm.mod3),data=working,x=T,y=T),working$iso3cc)$var

lm.mod4<-step(lm.mod3,steps=my.steps)
lm.mod4$se<-robcov(ols(formula(lm.mod4),data=working,x=T,y=T),working$iso3cc)$var

## NOW WITH VUL
lm.mod5<-update(lm.mod1,vul~.)
lm.mod5$se<-robcov(ols(formula(lm.mod5),data=working,x=T,y=T),working$iso3cc)$var

lm.mod6<-step(lm.mod5,steps=my.steps)
lm.mod6$se<-robcov(ols(formula(lm.mod6),data=working,x=T,y=T),working$iso3cc)$var
#lm.mod1<-lm(defacto_t~Ourscore+rol+runvetops+runpolarization+log(runpop)+log(rungdpch)+log(age)+HallGingerichCoordination+Laver.Hunt+log(share)+Turnover+Staff,data=working) #,x=T,y=T)

copyme<-apsrtable(lm.mod1,lm.mod2,lm.mod3,lm.mod4,lm.mod5,lm.mod6,
	se="robust",
	stars="default",
	model.names=c("De facto.1","De facto. 2","TOR. 1","TOR. 2","VUL. 1","VUL. 2"),
	coef.names=c("Intercept","De jure","Rule of law","Veto players","Polarization","log(Population)",
	"log(GDP per cap.)",
	"log(Age)","Coordination",
	"Salience",
	"Turnover","Staff"),
	Sweave=TRUE
)
cat(as.character(copyme),file="ols_models.tex",append=FALSE)



### 

library(Zelig)
z<-zelig(formula(lm.mod2),data=working,model="ls")
first_differences<-NULL

## Population
x.out.low<-setx(z,runpop=mean(working$runpop)-0.5*sd(working$runpop,na.rm=T))
x.out.high<-setx(z,runpop=mean(working$runpop)+0.5*sd(working$runpop,na.rm=T))
s.out<-sim(z,x.out.low,x.out.high)
first_differences<-rbind(summary(s.out)$qi.stats$fd,first_differences)  ## 0.17
## Rule of law
x.out.low<-setx(z,rol=mean(working$rol)-0.5*sd(working$rol,na.rm=T))
x.out.high<-setx(z,rol=mean(working$rol)+0.5*sd(working$rol,na.rm=T))
s.out<-sim(z,x.out.low,x.out.high)
first_differences<-rbind(summary(s.out)$qi.stats$fd,first_differences)  ## 0.17
## Veto players
x.out.low<-setx(z,runvetops=mean(working$runvetops)-0.5*sd(working$runvetops,na.rm=T))
x.out.high<-setx(z,runvetops=mean(working$runvetops)+0.5*sd(working$runvetops,na.rm=T))
s.out<-sim(z,x.out.low,x.out.high)
first_differences<-rbind(summary(s.out)$qi.stats$fd,first_differences)  ## 0.17
## De jure
x.out.low<-setx(z,Ourscore=mean(working$Ourscore)-0.5*sd(working$Ourscore,na.rm=T))
x.out.high<-setx(z,Ourscore=mean(working$Ourscore)+0.5*sd(working$Ourscore,na.rm=T))
s.out<-sim(z,x.out.low,x.out.high)
first_differences<-rbind(summary(s.out)$qi.stats$fd,first_differences)  ## 0.17

rownames(first_differences)<-c("Formal ind.","Veto players","Rule of law","Population")
first_differences<-as.data.frame(first_differences)
names(first_differences)[1]<-"Change"
names(first_differences)[3]<-"Lower"
names(first_differences)[4]<-"Upper"
first_differences$Variable<-rownames(first_differences)
first_differences$Variable<-factor(first_differences$Variable,
	levels=first_differences$Variable[order(first_differences$Change)],
	ordered=TRUE)
pdf(file="first_differences.pdf",height=2,width=5)
se<-ggplot(first_differences, aes(Variable, Change, 
  ymin = Lower, ymax=Upper),fontsize=8)
se + geom_pointrange() + coord_flip()
dev.off()

